* Gary LaMotte's Mathematica solutions to exercises from the 2012 Gravitational Wave Astronomy Summer School

* DAY1     PART4. FILTERING*

*PART a.*

samprate = 16384;

nsec=1.0;

ndata=nsec*samprate;

time = Table[(r-1.0)/ndata  , {r,ndata}];

sines = 5.*Sin[2.*Pi*74.7*time] + 2.*Sin[2.*Pi*820.*time];

noise = Table[RandomVariate[NormalDistribution[0.,1.4],ndata]];

ts = noise + sines;

fixed = Drop[ts,-1];

shiftedbyone = Insert[fixed,fixed[[1]],1];

*BELOW IS THE SIMPLE  FIRST-DIFFERENCE  FILTER*

subfilt = ts - shiftedbyone;

st = Transpose[{time, subfilt}];

ListLinePlot[st, PlotLabel -> Filter Difference, AxesLabel ->  {Time (Secs), Amplitude}]

ListLinePlot[st, PlotLabel -> Filter Difference, AxesLabel ->  {Time (Secs), Amplitude}, PlotRange-> {{0,0.05}, {-7,7}}]

ftsines = Fourier[subfilt];

ListLogPlot[(Abs[(ftsines)])^2,PlotRange->{{0,16384},{5,300000}},Joined->True, PlotLabel-> Periodogram, AxesLabel-> {Frequency, Amp^2}]

ListLogPlot[(Abs[(ftsines)])^2 ,PlotRange -> {{0,1000},{.01, 300000}},Joined->True, PlotLabel-> Periodogram, AxesLabel-> {Frequency, Amp^2}]

*PART c.*

*HAMMING FILTER*

hamfilt = (0.54 - 0.46(Cos[2*Pi*time]))*subfilt;

graph = Transpose[{time, hamfilt}];

ListLinePlot[graph, PlotLabel ->  Filter Difference Hamming, AxesLabel ->  {Time (Secs), Amplitude}, PlotRange-> {{0,1.}, {-10,10}}]

ftfilt = Fourier[hamfilt];

nyqpart = ftfilt[[1;;8192]];

ListLogPlot[(Abs[nyqpart])^2,PlotRange -> {{0,8192},{.01, 100000}}, Joined->True, PlotLabel-> NonSegmented Hamming,
AxesLabel-> {Frequency, Amp^2}]

ListLogPlot[(Abs[nyqpart])^2,PlotRange -> {{0,1024},{.01, 100000}}, Joined->True, PlotLabel-> NonSegmented Hamming, AxesLabel-> {Frequency, Amp^2}]

ListLogPlot[(Abs[nyqpart])^2,PlotRange -> {{0,130},{.00001, 100}}, Joined->True, PlotLabel-> NonSegmented Hamming, AxesLabel-> {Frequency, Amp^2}]

*BELOW: POWER SPECTRA WITH WELCH METHOD
PURPOSE:  DECREASE NOISE WITH UNFORTUNATE SPECTRAL LEAKAGE*

*Set up 8PARTITIONS and new time length*

partitions = 8;

segments = Partition[subfilt,     Length[subfilt]/partitions,      Length[subfilt]/(2*partitions)];

nsec2 = nsec / partitions

ndata2 = (samprate * nsec) / partitions

newtime = Table[(r-1.0)*(nsec2/ndata2)  , {r,ndata2}];

*DO HAMMING FILTERING OF EACH SEGMENT BELOW*

filtsegs = Transpose[(0.54 - 0.46(Cos[2*Pi* newtime/nsec2]))*Transpose[segments]];

*Do Hann filtering instead if desired: filtsegs = Transpose[Transpose[segments] * (1-Cos[2*Pi*newtime/nsec2])];*

f = filtsegs[[13]];

graph2 = Transpose[{newtime, f}];

ListLinePlot[graph2, PlotLabel -> Filter Difference Hamming Segmented, AxesLabel ->  {Time (Secs), Amplitude}, PlotRange->{{0,.125},{-10,10}}]

ftfiltsegs = (Abs[Fourier[filtsegs]])^2;

added= Total[ftfiltsegs];

averaged = added/((partitions*2) - 1);

index = Table[(r-1.0)*partitions  , {r,Length[averaged]}];

normalized = Transpose[{index, averaged}];

ListLogPlot[normalized,PlotRange -> {{0,8192},{.01, 100000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch, AxesLabel-> {Freq, Amp^2}]

ListLogPlot[normalized,PlotRange -> {{0,1024},{.01, 100000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch, AxesLabel-> {Freq, Amp^2}]

ListLogPlot[normalized,PlotRange -> {{0,130},{.00001, 100}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch, AxesLabel-> {Freq, Amp^2}]

*DAY 1.    4.FILTERING   (d)

*BELOW:  SIMULATED BUTTERWORTH FILTER BY USING SQUARE BANDPASS FILTER 65 TO 85Hz*

*SUCCESSFUL SQUARE WINDOW FILTER DESIGN IN FREQUENCY DOMAIN,INVERSE FOURIER TRANSFORMED TO THE TIME DOMAIN (SINC FUNCTION)AND CONVOLVED  WITH TIME SERIES*

*SQUARE FILTER PART*

start = {1,1};

pass = PadRight[start, 21, 1];

rtside = PadRight[pass, 16320];

sqfilt = PadLeft[rtside,16384];

*END SQUARE*

sinc = InverseFourier[sqfilt];

realsinc = (1./64.) * Re[sinc];

sincseries = Transpose[{time, realsinc}];

ListLinePlot[sincseries]

convts = ListConvolve[realsinc,ts,1];

sincts = Transpose[{time, convts}];

ListLinePlot[sincts, PlotLabel -> Filter BandPass, AxesLabel ->  {Time (Secs), Amplitude}, PlotRange-> {{0,1.}, {-6,6}}]

ListLinePlot[sincts, PlotLabel -> Filter BandPass, AxesLabel ->  {Time (Secs), Amplitude}, PlotRange-> {{0,.3}, {-6,6}}]

ftconvts = Fourier[convts];

ListLogPlot[(Abs[(ftconvts)])^2 ,PlotRange -> {{0,8192},{2, 300000}},Joined->True, PlotLabel-> Periodogram, AxesLabel-> {Frequency, Amp^2}]

ListLogPlot[(Abs[(ftconvts)])^2 ,PlotRange -> {{0,150},{2, 300000}},Joined->True, PlotLabel-> Periodogram, AxesLabel-> {Frequency, Amp^2}]

*DAY 1   PART 5. DOWN-SAMPLING, ANTI-ALIAS  FILTER*

*Part a.*

sixteenth = Partition[ts,1,16];

downsamp = Flatten[sixteenth];

Length[downsamp]

StandardDeviation[downsamp]

StandardDeviation[ts]

*Standard deviations are very similar btwn whole timeseries and 1/16th downsampled timesseries.*

*Part  b.*

*LOW PASS FILTER PART (512 Hz CORNER FREQ).  FILTER MADE FROM A SIMPLE SQUARE WINDOW FILTER IN FREQ DOMAIN, INVERSE FOURIER TRANSFORMED TO A SINC FILTER IN THE TIME DOMAIN, THEN CONVOLVED WITH TIME SERIES.*

start = {1,1};

pass = PadRight[start, 512, 1];

lowpass = PadRight[pass, 16384];

sinc = InverseFourier[lowpass];

realsinc = (1./64.) * Re[sinc];

lowpassfiltered = ListConvolve[realsinc,ts,1];
*THIS APPLIES THE LOW PASS FILTER(512Hz CORNER FREQ)TO THE TIME SERIES*

*NOW SELECT EVERY 16TH SAMPLE*

sixteenth = Partition[lowpassfiltered,1,16];

downsampfiltered = Flatten[sixteenth];

StandardDeviation[downsampfiltered]

StandardDeviation[ts]

*THE TIME SERIES WITH THE ANTI-ALIAS FILTER HAS A LOWER STD DEVIATION THAN THE TIME SERIES WHERE ALIASING IS PERMITTED (DUE TO THE 820 Hz SINE WAVE) AND A LOWER STD DEV THAN THE ORIGINAL TIME SERIES.*